Effect of Surface Andreev Bound States on the Bean-Livingston Barrier in d-Wave 

Superconductors 
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We study the influence of surface Andreev bound states in d-wave superconductors on the Bean- 
Livingston surface barrier for entry of a vortex line into a strongly type-II superconductor. Starting 
from Eilenberger theory we derive a generalization of London theory to incorporate the anomalous 
surface currents arising from the Andreev bound states. This allows us to find an analytical expres- 
sion for the modiflcation of the Bean-Livingston barrier in terms of a single parameter describing 
the influence of the Andreev bound states. We find that the field of first vortex entry is significantly 
enhanced. Also, the depinning field for vortices near the surface is renormalized. Both effects are 
temperature dependent and depend on the orientation of the surface relative to the d-wave gap. 
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Formation and dynamics of vortices is of interest in 
such diverse fields as cosmology, liquid crystals, and su- 
perfluids [!]■ It is well known that fermionic states ex- 
isting inside a vortex can significantly alter the structure 
of a vortex line [2^. In unconventional superconductors 
fermionic states, so-called Andreev bound states, can also 
appear at the surface. In this work we show that such 
surface states strongly influence the entrance and pin- 
ning properties of a vortex. Entry of a vortex into a 
type-II superconductor is hindered by the so-called Bean- 
Livingston (BL) surface barrier This barrier arises 
due to a competition of two forces acting on the vortex 
line: a force coming from the external magnetic field driv- 
ing the vortex into the superconductor and a force of the 
vortex' mirror image attracting it towards the outside. 
As a result, the penetration field Bg of first vortex entry 
is typically much larger than the lower critical field Bd 
and becomes of the order of the thermodynamic critical 
field Be i- 

At the surface of a d-wave superconductor Andreev 
bound states exist, depending on the relative orientation 
of the d-wave gap function to the surface @, 0, 0) Q • It 
has been shown that these states split in the presence of 
a supercurrent running along the surface, generating a 
quasi-particle current directed opposite to the supercur- 
rent [9]. At low temperatures this anomalous Meissner 
current even over-compensates the supercurrent, leading 
to a reversal of the current flow at the surface [9, 10]. 
Even though the Andreev bound states only exist within 
a coherence length ^ from the surface, they are influenced 
by a vortex line already, when the vortex is a penetra- 
tion depth A 3> ^ away from the surface, because the 
supercurrent field around the vortex is long-ranged 11 1 . 
Therefore, a modification of the BL barrier has to be 
expected due to the change of the free energy of the An- 
dreev bound states as a function of the position of the 



vortex creating an additional force. Here, we study this 
modification of the BL barrier and show that it results 
in a sizeable increase of the penetration field and the de- 
pinning field. 

We consider a superconducting half-space in the re- 
gion a; > 0. Both the cylindrical Fermi surface of the 
d-wave superconductor and the external magnetic field 
Bq shall be aligned parallel to the z-axis. In the follow- 
ing the single Abrikosov vortex is situated at the position 
r„ — (xy, 0) on the a:;-axis, and the specular surface shall 
be given by the y-axis. Since we wish to explore An- 
dreev bound states in an inhomogenous superconductor, 
Ginzburg-Landau theory is not sufhcient and as a start- 
ing point the minimum theory necessary is Eilenberger 
theory 12, [l3, fl^ - In particular, we use the Riccati- 
parametrization |l5l. [lij . where local equilibrium prop- 
erties of the d-wave superconductor are determined in 
terms of two scalar coherence functions a and b. To each 
unit vector k parametrizing the Fermi surface of the su- 
perconductor corresponds a quasiparticle trajectory in 
real space which is parallel to the Fermi velocity vj;-(k). 
Along this trajectory, and for given Matsubara frequency 
e„ — {2n+ l)iTkBT, two decoupled differential equations 
of the Riccati type have to be integrated 

(;iVF,V)a+(2£„-2ie/c(vF,A) + A*a)a- A = C(l) 
(avF,V)6-(2e„-2ie/c(vF,A)+A6)6 + A* = 0. 

Here, the arguments have been omitted for brevity, and in 
the following we use the specific gauge, in which the gap 
function A is real and a phase gradient is absorbed into 
the magnetic vector potential A. Starting from known 
bulk values, solutions for a (6) can be calculated in a 
stable way numerically, if the integration of the corre- 
sponding Riccati equation is performed parallel (antipar- 
allel) to the trajectory direction. Both the gap function 
A(r, k) — AooX(k)?/'(r) and the magnetic vector poten- 
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tial A(r) can be regarded as spatially dependent input 
parameters for the Riccati equations ([T]). In general, how- 
ever, their correct selfconsistent form has to be deter- 
mined from two additional conditions. One of them is 
given by the gap equation 

/ 2 \ 

Ao.i;{r) ^2nNoVkBT ^ / ^(k)— ^(r, k, e„) V 

e„>0 ^ ' 

(2) 

where A'o is the density of states in the normal phase, V 
the interaction strength, lOc the cut-off energy and (...) 
denotes an average over the Fermi surface of the super- 
conductor. Furthermore, the momentum dependence of 
the d-wave gap is implemented by x(k) = cos2(/). The 
other condition regards electrodynamic selfconsistency. 
It can be ensured via the quasiclassical expression for 
the current density 

i{r) = -AmeNoksT /v^(k)i_^(r, k, e„)\ (3) 

together with the Maxwell equation 

V X V X A(r) ^ — j(r). (4) 

c 

The system of equations is self-contained. In prin- 

ciple it can be used to find a fully selfconsistent solution, 
which is very time-consuming, however. 

In the following we present a way to derive a very 
good approximation to the fully selfconsistent solution of 
Eqs. We concentrate on a d-wave superconduc- 

tor with K ~ X/(, ^ 1. Furthermore the superconductor 
is considered to be clean, but not superclean. In other 
words, the effective mean free path vpT along a quasipar- 
ticle trajectory may be larger than the coherence length 
f — hvp/Aao, but it should also hold A ^ vpT. As a 
starting point we then use the selfconsistent solution of 
Eqs. ((I])-© in the absence of any magnetic vector poten- 
tials and currents. This solution can be found numeri- 
cally comparatively easy since the problem is translation- 
ally invariant along the boundary. This solution already 
contains both the existence of low-energy Andreev bound 
states and the local suppression of the gap in the vicin- 
ity of the boundary. In a next step, we start perturbing 
this known selfconsistent solution by a magnetic vector 
potential. Due to the assumptions stated above, we can 
restrict ourselves to perturbations, which typically vary 
on the length scale A and are effectively homogeneous 
along the mean free path vpT. In that case, the Riccati 
equations ([T]) can be evaluated for a constant Doppler 
shift in energy, ^{vf,A), along the quasiparticle trajec- 
tory. Let a denote the angle between the surface normal 
and a maximum gap direction of the d-wave (cf. inset 
of Fig. [T]). For the two most symmetric cases a = 
and a = 7r/4 respectively, the expansion of Eq. ([3]) to 
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FIG. 1: (Color online) Dimensionless function /, determining 
the local current response according to Eq. ([5}, as a function 
of the distance x to the boundary of a d-wave superconductor. 
Curves are shown for different temperatures and orientations. 
London electrodynamics, i.e. / = 0, holds for the d-wave 
orientation a — and towards the bulk. For a = n/A the 
current response in the vicinity of the boundary is modified 
due to surface Andreev bound states. This efi'ect is strongly 
enhanced for lower temperatures. 

first order of the perturbing vector potential A yields an 
equation of the form 

iir)-~^^[l-fir)]Air), (5) 

where the dimensionless function / on the righthand side 
is independent of A and can be determined numerically 
from the known fully selfconsistent solution for A = 
mentioned above. 

Eq. ^ presents a generalization of London theory in- 
cluding surface effects. Physically, the influence of both 
Andreev bound states and a gap suppression at the sur- 
face on the electrodynamics of the superconductor are 
completely represented by the function /. Due to transla- 
tional symmetry / only depends on the distance x to the 
boundary. In Fig. [T]we show our results for f{x). For the 
d-wave orientation a — neither Andreev bound states 
nor a suppression of the gap exist at the boundary and / 
vanishes identically. In this case London electrodynam- 
ics is valid everywhere in the superconducting half-space. 
For the orientation a = 7r/4 however, there are significant 
changes of the current response close to the boundary, 
since the transport is dominated by anomalous Meissner 
currents of the zero-energy Andreev bound states. The 
net current response 1 — / may even become negative, 
and this effect is enhanced for lower temperatures since 
the spectral weight of Andreev bound states strongly in- 
creases [1, A few coherence lengths into the bulk 
London electrodynamics is recovered. An important as- 
pect of the resulting Eq. ([5]) is, that the thermodynamic 
and electrodynamic selfconsistency conditions are effec- 
tively decoupled. As mentioned above, the function / can 
be obtained using the gap equation ^ only, without any 
electrodynamics considered. On the other hand, once 
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/ is provided, the magnetic vector potential A of the 
fully selfconsistent problem may simply be found from 
Eq. ^ and the modified London equation together 
with Neumann boundary conditions. 

In the following, we show how to find the resulting 
magnetic field and derive the modified surface barrier. 
For this purpose it is useful to employ the subscript no- 
tation B = Bi + Bjj, where the subscript L refers to 
the well-known London solution of the problem, whereas 
R denotes the unknown small 'response' term gener- 
ated by the function / in Eq. (O due to the presence 
of Andreev bound states. The London solution itself 
is given by Bl = Bf + Bl with Bf{r) = Boe~='/^z 



and Bl 



(r) 



Ko (^) - Ko (^) 



z where 



$0 is the flux quantum, Kq is a modified Bessel func- 
tion and the star operator shall define a mirror opera- 
tion with respect to the boundary, i.e. r* = {—x,y). 
Via the superscript notation a field is uniquely divided 
into two parts. The superscript W refers to the 'wall' 
part, which is translationally invariant along the bound- 
ary and is due to the screened external magnetic field. 
The 'vortex' part denoted by V incorporates the fields 
due to vortex and its mirror antivortex. The correspond- 
ing magnetic vector potential of the London solution is 
found from Al = -A^V x B^. Then, Eqs. (g]) and ^ 
may be combined to 



I 

A2 



Bi?(r) 



1 

A2 



V X [/(x)A(r)] . (6) 



Inversion of the differential operator leads to an integra- 
tion over the superconducting area S that can be written 
in the form 

Bi^W = ^ ^ d\'f[x') [V'Gs(r, r')] x A(r'), (7) 

where the appropriate Green's function Gs of the prob- 
lem has been involved, which is given by Gs (r,r') = 



Kq ^^^"x^^ ^ Ko {j^^ — '■^ • particular, the un 
known magnetic field B^,; according to Eq. ([7]) fulfils 
Eq. ([6]), is divergence- free and vanishes everywhere at 
the boundary. The integrand in Eq. ^ is weighted by 
the function /, which is effectively non-zero only up to a 
finite distance If <^ X from the boundary, where the stan- 
dard London theory is not valid (cf. Fig. [T]). Further- 
more, also the total deviation from standard London the- 
ory is small, i.e. J dx f X. Thus, concentrating on the 
main contribution to B/j, the source term A in the inte- 
grand can be replaced by the London solution A/,. More- 
over it is sufficient to keep the dominating term only: 
Sfl(r) = ^ Jsd^r'fix') [d^,Gsir,r')] AlJv'). This for- 
mula allows to calculate, how the magnetic field deviates 
from the London solution due to the presence of Andreev 
bound states. 

Let us derive the modified BL barrier now. The corre- 
sponding free energy density g of our system is given by 



9 = 9o + ^[j^{l~ /)A2 -f B2 - 2Bo • B] . Here, go de- 
notes the value in the absence of magnetic fields and cur- 
rents. A rather long but straightforward analysis yields 
the resulting total free energy G (per unit length of the 
z-direction) as a function of the vortex distance Xy : 



_$o_ 



Koi2x,/X) 



^0 Js 



(8) 



Here, an arbitrary additive constant has been chosen such 
that G{xy oo) = 0. Note, that the first two terms are 
the original BL barrier (cf. The latter terms are 

modifications which only appear when f ^ and thus 
Bj^ ^ 0, correspondingly. After some algebra, we find 



""/^ r dx'fix') (l + e-2-7A^ 



Bl{v.) = dx'fix')K{xy,x'), 

and B^{rv) — d^r_B]((r) where the abbreviation 

K{x,,x')^Ki (2^^'^+2Ki{^)+Ki (2^^) has 

been used and we restricted ourselves to Xy >lf, i.e. the 
vortex is not directly situated in the small surface re- 
gion with / ^ 0. Since the function / is nonzero only 
in a small surface region, the above expressions can be 
further simplified, if one replaces f{x') by a delta func- 
tion CfS{x') with Cf = dx' f{x'). The quantity Cf has 
the dimension of a length and remains the only parame- 
ter describing the influence of the Andreev bound states. 
Then, the final result for the normalized free energy of 
Eq. dH) is 



G{xy) = 2Bo(I + c/)e- 



[Ko{2xy) + 2cfKi{2xy)] 



(?) 

Here, the normalized G is defined as G/($Q/I67r^A^), Bq 
denotes the applied magnetic field in units of $o/27rA^, 
and we abbreviate Xy = Xy/X and c/ = Cf/X. 

Eq. ([9]) is our generalization of the surface barrier. 
For Cf = it reduces to the result of Bean and Liv- 
ingston [3, In the presence of Andreev bound states, 
however, there are two modifications due to c/ ^ 0. The 
first modification is an enhancement of the repulsive part 
of the potential, effectively increasing the external mag- 
netic field amplitude. Secondly, also the attractive part 
of the potential coming from the antivortex is strength- 
ened by an additional term. Note, that c/ <C f. Even 
the maximum value of Fig. [T]occuring for the orientation 
a — 7r/4 at the temperature T — Q.lTc is only c/ w 3.7If , 
i.e. Cf ~ 3.71k~^. Nevertheless particularly the second 
term becomes important, because it contains a Bessel 
function, which diverges for small vortex distances Xy. 

It is important to realize, how drastic some proper- 
ties may change just because of this rather small looking 
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FIG. 2: (Color online) Relative change of the external mag- 
netic field B+ needed to release a vortex from a pinning cen- 
ter towards the bulk. Results are shown as a function of the 
normalized pinning force Fp for different positions Xp of the 
pinning center. Obviously, unpinning the vortex in the pres- 
ence of surface bound states requires a significantly higher 
external magnetic field. 



FIG. 3: (Color online) Temperature dependence of the pa- 
rameter for a = 7r/4. The inset shows the expected 
temperature dependence of the field of first vortex entry Bs 
according to Eq. (111(1 (solid line) compared with the BL result 
(dashed line). 



modification. To give an example we consider a vortex 
pinned in a pinning center at the distance Xp from the 
boundary. The maximum pinning force of that center, 
measured in units of ^q/Stt^A^, is characterized by the 
parameter Fp. If the applied magnetic field then exceeds 
some specific value 5+, which can be determined from 
Eq. ([HI), the vortex gets unpinned and proceeds towards 
the bulk. The relative change for this depinning field is 
given by 



B+{cf) 
B+{0) 



1 



_ /r,;; Ko(2x^) CJ_ 



1 



(10) 



For the parameters k — 200 and c/ = 3.71^ from above, 
some results are shown in Fig. [2l Clearly, the external 
magnetic field needed to unpin a vortex close to a bound- 
ary which exhibits Andreev bound states is significantly 
enhanced in the most relevant cases. 

Analogous to the classical BL work [sl , we can roughly 
estimate the field of first free vortex entry Bg by choosing 
Fp = and a distance Xy ^ ^ from the boundary. Then, 
we find 



Bs{0) 



1 



21 



(11) 



where we have assumed c/ <C 1, but c//^ ~ 1. Here, 
Bs{0) is the BL value. This result means that the field 
of first vortex entry is significantly enhanced by the pres- 
ence of Andreev bound states. Compared to the classical 
BL value it may increase by a factor 4-5 for low temper- 
atures and a = 7r/4. In Fig. [3] we show the temperature 
dependence of c//^. The corresponding characteristic 
change of Bs as a function of temperature, which is indi- 
cated by the inset of Fig. [31 may facilitate experimental 
observation of this effect. 



We also expect the modification of the BL barrier to 
have a strong influence on the hysteresis of the split of 
the zero bias conductance peak as observed in tunneling 
However, a detailed calculation of 
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experiments 

the hysteresis curve goes beyond the scope of the present 
work. 
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